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Phase  separation  within  water  containing  polymer  membranes  is  studied  by  dissipative  particle  dy¬ 
namics  (DPD).  The  polymers  are  composed  of  hydrophobic  backbone  A  beads  to  which  ([A-C]  or  [A-A-C- 
C])  side  chains  are  attached  with  a  hydrophilic  (C)  end  moiety.  Water  diffusion  through  the  water 
containing  pores  is  modeled  by  Monte  Carlo  (MC)  tracer  diffusion.  Several  polymeric  architectures  are 
considered  which  differ  in  the  way  side  chains  are  attached  along  the  polymer  backbones.  In  total  120 
pore  morphologies  are  stored  on  file  and  probed  by  MC  tracer  diffusion  calculations.  For  membranes  of 
the  same  ion  exchange  capacity  diffusion  is  highest  for  architectures  for  which  two  side  chains  are 
branching  off  from  the  backbone  sites  as  compared  to  architectures  for  which  only  one  side  chain  is 
branching  off.  This  conclusion  is  also  obtained  from  diffusion  constants  derived  from  mean  squared 
displacements  of  water  during  the  DPD  simulations.  The  intrinsic  high  polymer  mobility  in  conventional 
DPD  results  in  water  diffusivities  that  are  much  higher  than  those  obtained  from  the  MC  calculations. 
Assigning  higher  masses  to  the  polymer  beads  results  in  significant  decrease  in  water  bead  motion. 

©  2014  Published  by  Elsevier  B.V. 


1.  Introduction 

Nafion®  is  applied  in  water  purification  by  reverse  osmosis  1  , 
dialysis  and  as  a  proton  exchange  membrane  (PEM)  in  polymer 
electrolyte  fuel  cells  (PEFC).  Within  a  PEM  Nafion  separates  the 
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cathode  from  the  anode  and  simultaneously  acts  as  a  conductor  for 
protons.  The  chemical  formula  of  Nafion  is  given  in  Fig.  1(a).  Nafion 
is  a  per-fluorinated  sulfonic  acid  (PFSA)  ionomer  that  has  side 
chains  attached  to  its  hydrophobic  Teflon  (— CF2— )  backbone.  Since 
the  side  chains  contain  pendant  hydrophilic  SO3H  groups  the 
membrane  takes  up  water  when  exposed  to  humid  environments. 
The  water  uptake,  expressed  as  the  number  of  water  molecules  A 
per  sulfonic  site,  increases  with  relative  humidity  [2-9].  Phase 
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Fig.  1.  (a)  Repeat  unit  of  Nafion®(EW  =  1143)  and  DPD  parameterization  as  proposed 
in  Ref.  42.  A':  CF2-CF2-CF2-CF2i  B' :  0-CF2-CF(CF3)-0,  C:  CF2-CF2-S03H.  (b)  Repeat 
unit  of  the  polymer  types.  Side  chains  are  branching  off  from  every  4th  (I)  or  8th  (III, 
IV)  A  bead  along  the  backbone.  For  type  II  side  chain  branching  occurs  at  every  7th  and 
8th  backbone  A  bead.  The  number  of  bonds  or  topological  distance  Dtopoi  i  between 
nearest  C  beads  are  given  for  types  II,  III,  and  IV. 


separation  results  in  a  network  with  water  containing  pores  a  few 
nm  in  diameter  surrounded  by  polymer  phase  [7-9  .  Because  the 
pores  provide  the  pathways  for  water  and  proton  diffusion,  the 
degree  of  pore  connectivity  is  important  to  obtain  high  proton 
conductivities.  At  low  humidity  the  proton  conductivity  decreases 
rapidly  2-5,7,10,11],  which  requires  a  fuel  cell  design  strategy  that 
provides  sufficient  hydration  of  the  PFSA  membranes  during 
operating  conditions. 

Much  effort  is  spent  on  the  development  of  new  polymers  with 
high  proton  conductivity.  These  alternatives  have  acidic  sites 
located  within  the  backbones  of  the  block  polymers  or  within  the 
side  chains  of  the  grafted  polymers,  such  as  block  polymer  type 
poly-benzimidazole  [12],  sulfonated  poly-ether-ether-ketone 
[13,14  ,  sulfonated  styrene/ethylene-butylene/styrene  [10]  sulfo¬ 
nated  poly  aryl  ether  ketone  [15]  and  bi-phenyl  sulfone  [16  ,  sul¬ 
fonated  polyimide  [17,18],  and  the  grafted  [bis]  (perfluoroalkyl) 
sulfonyl  imide  perfluorinated  ionomer  [19  .As  for  the  PFSA  mem¬ 
branes  the  presence  of  hydrophilic  acidic  sites  cause  membranes 
composed  of  these  alternative  polymers  to  swell  and  phase  sepa¬ 
rate  with  water  uptake. 

The  amount  of  water  uptake  and  pore  connectivity  depends 
on  many  parameters  such  as  molecular  architecture,  side 
chain  distribution  and  side  chain  length  (grafted  polymers), 
acidic  site  density,  backbone  and/or  side  chain  stiffness,  and 
molecular  weight.  Molecular  dynamics  (MD)  18-31],  coarse 
grained  MD  (CGMD)  [32,33]  and  meso-scale  approaches  such  as 
dissipative  particle  dynamics  (DPD)  [34-42]  can  provide  insight 
in  how  pore  size  and  pore  connectivity  are  affected  by 
these  parameters.  MD  has  been  applied  frequently  to  model 
phase  separation  within  PSFA  membranes  [18,20-31  .  These 
studies  predict  water  clustering  in  agreement  with  experiment. 
Water  and/or  proton  diffusion  constants  can  be  obtained  from 
the  mean  squared  displacements  (MSD)  of  these  species.  The 
period  over  which  sampling  is  performed  in  MD  is  typically  a  few 
ns  which  is  too  small  to  probe  a  representative  pore  network 
within  real  membranes.  That  would  require  a  meso  scale  simu¬ 
lation  in  which  the  system  size  is  significantly  larger  than  the 
characteristic  distance  (eg.  inter  cluster  distance)  of  the  pore 
network.  Since  within  PFSA  membranes  the  characteristic  dis¬ 
tance  is  around  5  nm  [8,9,11]  obstructions  in  the  diffusion 
pathways  can  be  probed  only  for  system  volumes  of  O(104  nm3) 
or  more. 


For  system  size  of  -60  nm3  Devanathan  et  al.  27]  studied  with 
MD  water  and  proton  diffusion  through  Nafionll43  (1143  is  the 
equivalent  weight  (EW)  of  1143  g  polymer/mol  SO3H).  A  significant 
increase  in  diffusion  occurred  near  A  =  5.  For  three  times  larger 
system  size  Cui  et  al.  [28]  also  observed  increased  water  and  proton 
diffusion  within  Nafion  with  water  uptake,  in  accordance  with 
experiment.  For  a  volume  of  -400  nm3  Karo  et  al.29  [29]  simulated 
water  diffusion  within  Dow  (whose  side  chains  are  composed  of 
O-CF2-CF2-SO3H)  and  Nafion  at  A  =  15.  They  noted  that  the  pore 
morphology  depends  strongly  on  system  size  and  questioned 
conclusions  [30]  obtained  from  smaller  system  size.  An  MD  study 
by  Knox  and  Voth31  of  Nafion  that  contained  -20  vol%  water 
involved  2  x  106  atoms  (-27000  nm3).  Since  the  evolvement  from  a 
random  configuration  to  an  equilibrium  structure  is  impossible  for 
such  large  system,  they  pre-assumed  5  different  morphologies  and 
a  random  morphology.  With  exception  of  the  random  morphology, 
each  of  them  evolved  such  that  the  characteristic  scattering  peak 
could  be  reproduced.  Therefore  solely  based  on  the  Bragg  reflection 
peak  position  in  scattering  experiments  no  unique  model  can  be 
assigned  to  the  specific  pore  morphology  within  Nafion. 

The  effect  of  side  chain  distribution  on  morphology  and  diffu¬ 
sion  was  studied  by  Jang  et  al.  [32]  by  CGMD.  They  found  less 
pronounced  phase  separation  and  narrower  pores  for  Nafion 
polymers  with  uniform  side  chain  attachment  as  compared  to 
polymers  with  all  side  chains  located  at  one  end  of  the  backbone. 
However,  water  and  proton  diffusion  was  not  drastically  affected  by 
these  differences  in  sequence  design.  Allahyarov  et  al.  [33] 
concluded  from  CGMD  simulations  that  for  PFSA  like  polymers 
with  similar  EW  those  with  the  longer  side  chains  form  larger 
clusters  and  higher  proton  conductivities.  This  is  consistent  with 
DPD  studies  [35,38]  obtained  for  larger  system  size. 

Dorenbos  et  al.  combined  DPD  with  Monte  Carlo  (MC)  tracer 
diffusion  calculations  to  predict  the  dependency  of  water  diffusion 
on  hydration  level  [34],  EW  [34,35],  side  chain  length  [34,35,38] 
and  Nafion  layer  thickness  (<10  nm)  on  substrates  of  various  hy- 
drophobicity  [36  .  The  main  idea  behind  the  DPD-MC  strategy  is 
twofold.  Firstly,  much  larger  equilibrium  morphologies  can  be  ob¬ 
tained  by  DPD  as  compared  to  MD  with  minimal  computational 
effort.  Secondly,  water  diffusion  can  be  simulated  separately  by 
performing  MC  tracer  diffusion  through  selected  (frozen  or  static) 
DPD  morphologies.  The  fast  evolvement  towards  equilibrium 
morphologies  is  caused  by  the  soft  repulsive  interaction  that  acts 
between  atomic  groups  that  are  coarse  grained  into  beads.  There¬ 
fore  topological  violations,  due  to  bond  crossings,  can  lift  up  en¬ 
tanglements  in  conventional  DPD,  which  explains  the  absence  of 
polymer  reptation  motions  (D-l  /N2,  with  N  number  of  polymeric 
backbone  beads)  [43].  Polymer  chains  diffuse  thus  much  faster  than 
in  an  MD  simulation. 

Due  to  the  polymer  bead  motion  the  pore  network  evolves 
while  water  is  moving  through  it.  This  results  in  higher  water  bead 
motion  through  the  dynamic  networks  as  compared  to  diffusion 
through  a  pore  network  in  which  polymer  motions  are  restricted  by 
entanglements.  This  might  in  a  DPD  simulation  obscure  the  de¬ 
pendency  of  solvent  (water)  diffusion  on  various  design  parame¬ 
ters.  Within  the  DPD-MC  scheme  this  problem  is  circumvented  by 
mapping  the  equilibrium-like  DPD  configurations  on  a  grid  and 
performing  MC  trajectory  calculations  through  the  mapped  pore 
networks.  In  a  recent  study  [41  it  was  verified  that  within  hydrated 
amphiphilic  block  polymers  diffusion  constants  derived  from  the 
water  bead  motions  were  indeed  higher  than  MC  derived  water 
diffusion  constants. 

Water  diffusion  constants  within  Nafion  derived  from  DPD-MC 
calculations,  in  which  water  diffusion  is  restricted  to  the  hydro¬ 
philic  phase  that  contain  both  water  and  sulfuric  acid  sites, 
approach  experimental  diffusion  values  [34  .  More  generally  it  was 


538 


G.  Dorenbos  /  Journal  of  Power  Sources  270  (2014)  536-546 


predicted  that  for  architectures  with  same  side  chain  lengths  an 
increase  in  EW  results  in  increased  water  cluster  size  and  distance 
between  them  but  with  slower  (water)  diffusion  [34,38].  Theoret¬ 
ical  studies  on  PFSA-like  polymers  in  which  also  side  chain  lengths 
were  varied  [34,35,38]  predicted  that  at  fixed  water  volume  frac¬ 
tion,  0W,  for  architectures  of  the  same  EW  those  that  contain  the 
longer  side  chains  exhibit  the  largest  water  clusters,  largest  inter 
cluster  spacing  and  fastest  diffusion.  For  the  longer  side  chain  ar¬ 
chitectures  also  lower  percolation  thresholds  are  expected  than 
those  containing  short  side  chains  [38  .  This  is  caused  by  different 
topological  distance  between  nearby  hydrophilic  sites  within  the 
architectures.  Architectures  with  larger  topological  distance  allow 
the  formation  of  larger  and  better  connected  pores,  higher  diffusion 
constants  at  fixed  </>w,  and  lower  percolation  thresholds.  A  modified 
kinetic  Monte  Carlo  based  DPD-MC  algorithm  was  developed  to 
study  gas  diffusion  through  membranes  in  which  gaseous  species 
have  access  to  both  polymer  and  water  phases.  This  permits  to 
predict  the  dependence  of  02  and  H2  gas  permeability  through 
hydrated  Nation  on  the  membranes'  EW  [37]  within  a  two  phase 
model. 

In  the  DPD  studies  addressed  above  side  chains  were  uniformly 
grafted  along  the  polymer  backbones.  Recently  it  was  predicted 
that  for  PFSA-like  ionomers  of  the  same  EW  a  non-uniform  side 
chain  distribution  results  in  better  connected  pores  with  increased 
diffusion  [39,40]  and  lower  percolation  thresholds  [40]  as 
compared  to  uniform  side  chain  attachments.  Also  the  agreement 
between  experimental  and  calculated  02  permeability  and  water 
diffusion  constants  through  hydrated  Nation  improved  for  the  more 
realistic  architectures  with  statistically  distributed  side  chains  [39]. 

The  purpose  here  is  to  predict  how  water  diffusion  through  the 
pore  networks  of  membranes  of  the  same  EW  might  be  affected  by 
location  of  the  side  chains  along  the  backbones.  Architectures  are 
considered  with  one  or  two  side  chains  branching  off  from  a  single 
backbone  bead,  and  side  chains  can  contain  one  or  two  hydrophilic 
sites  (Fig.  1(b)).  DPD  is  used  to  generate  pore  morphologies  at  a 
fixed  hydration  level  of  </>w  =  0.16.  Diffusion  through  the  pore 
networks  is  modeled  by  employing  Monte  Carlo  (MC)  tracer  cal¬ 
culations  through  selected  static  (or  frozen)  DPD  morphologies.  The 
diffusion  of  the  water  beads  through  the  dynamic  pore  networks  is 
also  studied  and  obtained  DPD  diffusion  constants  are  compared 
with  the  MC  derived  values. 

The  outline  of  this  paper  is  as  follows.  In  Computational  details 
the  polymer  architectures  are  presented  and  the  DPD  parameteri¬ 
zation  is  outlined.  In  Results  and  analysis  pore  morphologies  and 
relations  between  morphology,  chain  architecture  and  MC  derived 
diffusion  constants  are  presented.  It  is  found  that  differences  in 
diffusion  are  not  due  to  more  water  being  located  within  the 
percolating  clusters.  Diffusion  constants  derived  from  the  DPD 
water  bead  motions  are  compared  with  the  MC  derived  diffusion 
constants.  In  Discussion  the  limitations  of  the  used  computational 
approaches  are  addressed. 

2.  Computational  details 

2.1.  Polymer  chain  architecture 

In  DPD  several  atoms  or  molecules  are  grouped  into  a  single 
bead.  Four  types  of  polymers  are  considered.  A  single  repeat  unit  of 
each  type  is  shown  in  Fig.  1(b)  and  contains  hydrophobic  A  beads 
and  one  or  two  hydrophilic  C  beads.  For  type  I,  II,  and  III  the  side 
chains  are  represented  by  connected  [AC]  fragments,  whilst  a  side 
chain  for  type  IV  polymer  is  two  times  longer  and  represented  by 
[A2C2].  Side  chains  are  branching  off  from  every  4th  (I),  7th  and  8th 
(II),  and  8th  (III,  IV)  A  bead  along  the  backbone.  For  type  III  polymers 
two  side  chains  are  attached  to  a  single  A  bead.  The  DPD 


formulation  of  the  repeat  units  are:  A4[AC]  (I),  A7[AC]A1[AC]  (II), 
A8[AC][AC]  (III),  A8[A2C2]  (IV).  Six  architectures  are  considered,  see 
Table  1.  The  number  of  repeat  units  for  each  polymer  type  is  6 
(architectures  I,  Illb  and  IVb)  or  3  (architectures  II,  Ilia  and  IVa).  The 
number  of  side  chains  per  polymer  is  thus  NSide  =  3,  6,  or  12.  This 
choice  results  for  architectures  I,  II,  Ilia,  and  IVa  in  the  same  poly¬ 
mer  length,  ie.  identical  number  of  backbone  beads  (Nbackbone)  per 
polymer  with  Nbackbone  =  24.  The  side  chain  density,  pSide,  expressed 
as  the  ratio  of  the  number  of  side  chains  and  backbone  beads,  for 
type  I,  II,  and  III  are  the  same  with  pside  =  NsjdeNbackbone  =  1/4  (see 
Table  1 ).  For  type  IV  pSide  is  half  this  value.  Also  the  two  times  longer 
architecture  IVb  is  considered  with  the  same  number  of  side  chains 
as  type  I,  II,  and  Ilia.  For  architecture  IVb  Nbackbone  =  48,  for  this 
reason  also  architecture  Illb  is  included  with  Nbackbone  =  48  and 
Nside  =  12.  The  main  purpose  is  to  draw  trends  from  architectures 
for  which  the  lengths  are  the  same  (I,  II,  Ilia,  and  IVa).  Architectures 
Illb  and  IVb  are  included  to  verify  whether  the  same  trends  are 
observed  as  between  Ilia  and  IVa  for  increased  polymer  length. 

Water  is  represented  by  W  beads.  Within  the  DPD  framework 
followed  here  (see  next  paragraph)  the  volume  of  each  (A,  C  and  W) 
bead  is  the  same  and  equal  to  0.12  nm3  which  corresponds  to  the 
volume  occupied  by  4  water  molecules.  The  (dry)  acidic  site  den¬ 
sities  are  thus  2.31  mmol/cm3.  The  ion  exchange  capacity  (IEC) 
expressed  in  mmol/gr  can  be  obtained  when  the  mass  densities  of 
the  molecular  fragments  are  known.  For  instance,  for  a  polymer 
with  a  mass  density  of  1  gr/cm3  the  IEC  would  be  2.31  mmol/gr  (or 
EW  of  433  gr/mol). 


2.2.  Dissipative  particle  dynamics 


Hoogerbrugge  and  Koelman  44]  introduced  DPD  to  study  hy¬ 
drodynamics.  DPD  is  used  in  modeling  eg.  phase  separation  [45], 
membrane  rupture  [46  ,  vesicles  [47,48]  and  drug  delivery  49]. 
Here  the  framework  developed  by  Groot  and  Warren  50]  is  fol¬ 
lowed  which  can  be  considered  as  standard  and  therefore  the  most 
essential  details  are  given  here.  According  to  Newton's  equations  of 
motion  beads  evolve  as:  drz/dt  =  v;,  m[d\[\dt  =  fa  with  m*,  rz  and  v*, 
the  mass,  position  and  velocity  of  bead  i.  Interactions  between 
beads  i  and  j  are  given  by  conservative  Ffj ,  dissipative  F{/,  and 
random  Ffj  forces.  The  force /i  acting  on  bead  i  is  given  by  Eq.  (1). 

fi  EFV+FU+FU  U) 

j*i 

The  sum  is  over  all  particles  j  located  within  a  cutoff  distance  rc  that 
defines  the  length  scale.  Ffj  is  soft  repulsive  and  decreases  linearly 
with  distance  between  beads: 


(2a) 


Table  1 

Number  of  side  chains  (2nd  column),  backbone  beads  (3rd  column)  and  side  chain 
density  (4th  column)  for  each  architecture.  The  topological  distance  Dt0poi,i  (5th 
column)  and  Dtopoi,2,  (6th  column)  are  the  number  of  bonds  between  nearest  and 
next  nearest  C  beads  within  the  architectures,  respectively.  7th  column:  Total 
number  of  polymers,  iYpoi  in  simulation  box. 


Architecture 

Nside 

^backbone 

Pside 

f^topol,  1 

f^topol,  2 

l^pol 

I:  A4[AC]-6 

6 

24 

4_1 

8 

8 

968 

II:  A7[AC]Al[AC]-3 

6 

24 

4"1 

5 

11 

968 

Ilia:  A8[AC][AC]-3 

6 

24 

4"1 

4 

12 

968 

Illb:  A8[AC][AC]-6 

12 

48 

4"1 

4 

12 

484 

IVa:  A8[A2C2]-3 

3 

24 

8  1 

1 

14 

968 

IVb:  A8[A2C2]-6 

6 

48 

8  1 

1 

14 

484 
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(2b) 


In  Eq.  (2)  ?y  is  the  unit  vector  in  the  direction  of  bead  j  towards  bead 
i :  F?j  corresponds  to  thermal  noise  (Eq.  (3)),  and  fj/  is  proportional  to 
the  beads'  relative  velocity  (Eq.  (4)). 


F]j  =  «^(r#)fy(A0-a5(kiirr% 


if  =  -ycoD(rij)(rij-vij)?ij  (4) 

Vy  =  v,-Vj,  wD  and  coR  are  weight  functions,  introduces 
randomness  into  the  system  involving  a  randomly  fluctuating 
variable  with  Gaussian  statistics  with  zero  mean  and  unit  variance: 
<?y(t)  >  0,  and  <?y(t)?k/(t')  >  (Mj7  +  du5jk)5{t-t').  For  each  pair  of 
beads  there  is  an  independent  random  function.  coD  and  coR  are 
given  by  Eq.  (5). 


Noise  (a)  and  friction  (y)  parameters  are  related  as  a2  =  2y kBT 
[51  with  g  =  3,  y  =  4.5,  kB  the  Boltzmann  constant  and  T  the 
temperature.  All  three  forces  act  along  the  line  of  centers  and 
conserve  linear  and  angular  momentum.  Adjacent  beads  that 
belong  to  the  same  polymer  are  joined  by  a  spring  force  with  spring 
constant  C  =  50  and  equilibrium  bond  distance  Rq  =  0.85  rc: 


F§  =  -C(ry  -  R0)r y  (6) 

As  for  the  interaction  range  rc,  the  mass  of  all  beads  and  the  unit 
of  energy  (1<BT)  are  scaled  to  1.  The  unit  of  time  r  =  rc  (m//<Br)°-5  is 
also  equal  to  1.  The  temperature  is  kept  constant  near  1<BT  =  1.0  by 
solving  the  equations  of  motion  using  a  modified  Verlet  integration 
scheme  50]  with  empirical  factor  0.65  and  time  step  At  =  0.05.  The 
bead  density  p  is  3. 

The  repulsions  a y  are  listed  in  fable  2.  Those  between  incom¬ 
patible  (hydrophobic  vs  hydrophilic)  beads  are  highest  (A-W  and 
A-C).  Repulsions  between  similar  beads  are  set  at  104  which  re¬ 
produces  water  compressibility  [34,39,42  .  The  other  repulsions 
are  calculated  from  the  Flory-Huggins  %  parameters  ( Table  2)  that 
are  comparable  with  those  derived  from  calculated  binding  en¬ 
ergies  for  the  molecular  fragments  of  Nation  by  Yamamoto  and 
Hyodo  [42  .  They  parameterized  a  Nation  repeat  unit  as  A'4[B'C/] 
(Fig.  1(a))  with  hydrophobic  A'  (CF2-CF2-CF2-CF2)  and  B' 
(O-CF2-CF-CF3-O)  beads  and  hydrophilic  C  (CF2-CF2-SO3H) 
beads,  and  obtained  xa'w  =  5.79,  xb'w  =  4.9,  yew  =  -2.79, 
Xa'B'  =  0.02,  xa'c  =  3.11.  Here  polymers  are  constructed  of  hydro- 
phobic  A  and  hydrophilic  C  beads  (Fig.  1(b)).  The  %  parameters 
between  hydrophobic  and  hydrophilic  (C  and  W)  beads  are  chosen 
to  be  xac  =  Xaw  =  4.9.  This  value  is  equal  to  xb'w  =  4.9  [42]  and 
between  xaw  =  5.79  and  xa'c  =  3.11.  xcw  =  -2.6  is  close  to 
Xcw  =  -2.79  estimated  in  Ref.  [42]. 

The  y-values  in  Table  2  were  adapted  in  previous  works  38-41  ] 
and  resulted  in  well  phase  separated  morphologies.  They  are  also 
expected  to  result  in  phase  separated  morphologies  for  the  current 


Table  2 


DPD  repulsions  and  xy-parameters  (parenthesis). 


A 

C 

W 

A 

104 (0) 

C 

124.4  (4.9) 

104  (0) 

w 

124.4  (4.9) 

93.2  (-2.6) 

104  (0) 

II:  A7[AC]A1  [AC]-3 


IVa:  A8[A2C2]-3 


Fig.  2.  Snapshots  of  pore  morphologies  obtained  at  2  x  104At  for  four  membrane  types 
(C  beads:  yellow  (light  grey),  W  beads:  blue  (dark)).  (For  interpretation  of  the  refer¬ 
ences  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this 
article.) 


architectures.  For  the  sake  of  generality,  atomistic  calculations  to 
obtain  %  parameters  for  given  chemical  formulae  of  the  beads  were 
omitted.  Important  is  the  incompatibility  between  A  and  C  beads, 
with  (A)C  fragments  (un)favorably  interacting  with  water.  The  % 
parameters  in  Table  2  were  also  used  in  earlier  modeling  of  phase 
separation  within  amphiphilic  block  and  amphiphilic  grafted 
polymers.  Therefore  an  already  existing  database  is  extended  in  this 
work. 

The  DPD  repulsions  that  correspond  with  the  x-parameters  are 
obtained  from  Eq.  (7)  [41]. 

A  fly  =  (4.16±0.15)  x  Xij  (7) 

The  water  volume  fraction  0W  is  given  by  Eq.  (8),  with  NA,  Nw 
and  the  total  number  of  A,  W  and  C  beads. 

Nw 

=  j\jw  +  NA  +  Nc  ^ 

In  this  work  0W  =  0.16.  Since  each  W  bead  represents  4  water 
molecules,  the  water  content  expressed  in  terms  of  A  is  given  by 
Eq.  (9). 


Nw 

A  Nc  X  4 


For  0W  =  0.16  the  corresponding  A  value  is  4.6. 

The  system  is  cubic  with  cell  dimension  of  L  =  24.  One  DPD 
length  unit  corresponds  to  rc  =  (3  x  0.12)^3  nm3  =  0.71  nm.  Box 
edge  lengths  are  L  x  rc  ~17  nm  and  volumes  are  ~5  x  103  nm3.  The 
simulation  box  contains  41472  (=p  x  L3)  beads  (6624  W  beads  and 
Npoi  =  968  or  484  polymer  chains).  Periodic  boundaries  are  applied 
in  the  orthogonal  directions.  The  physical  time  to  which  r  corre¬ 
sponds  for  aww  =  104  is  -130  ps.  Morphologies  were  generated  up 
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to  4  x  104  At  and  correspond  with  a  duration  of  O  (0.2  ps).  with  diagonal  repulsions  of  aww  =  51  but  with  chi  parameters  being 

Figures  that  display  morphologies  were  generated  using  the  VMD-  the  same  as  here, 

visual  molecular  dynamics  package  [52  . 


3.  Results  and  analysis 


3.2.  MC  tracer  diffusion  through  static  pore  networks 


3.1.  Pore  morphology 

In  Fig.  2  snapshots  of  morphologies  obtained  at  t  =  2  x  104  At  are 
displayed.  Only  C  and  W  beads  are  shown.  For  all  morphologies,  C 
beads  are  associated  with  W  beads  and  form  a  hydrophilic  (C  +  W) 
pore  network.  For  architecture  I  the  pores  appear  to  be  most 
dispersed,  with  narrow  connections  between  the  clusters.  For  ar¬ 
chitecture  III  the  pores  are  most  well  developed  and  cylindrically 
shaped.  The  water  clusters  are  largest,  less  dispersed  and  most 
spherical  for  architecture  IV.  From  visual  inspection  of  the  mor¬ 
phologies  during  a  DPD  run  it  was  observed  that  the  pore  networks 
are  capable  to  slightly  change  their  morphology.  This  is  illustrated 
for  the  A8[AC][AC]-3  architecture  in  Fig.  3,  where  snapshots  taken 
at  intervals  of  4000At  (or  200 r)  are  displayed  (for  instance  by 
comparison  of  the  morphology  obtained  at  t  =  12000At  with  that  at 
t  =  32000At). 

The  extent  of  morphological  variation  during  the  DPD  runs  was 
also  examined  by  inspection  of  the  W-bead  pair  correlation  func¬ 
tions,  g(r).  Examples  of  the  time  evolution  of  g(r)  are  given  in 
Fig.  4(a)— (d).  The  average  distance  between  water  clusters,  Dci-ci, 
can  be  estimated  from  the  position  of  the  second  maximum  34,42] 
(located  between  4  nm  and  7  nm).  Dq-ci  is  plotted  against  DPD  time 
in  Fig.  5(a).  After  an  initial  increase  Dq-ci  tends  to  stabilize  beyond 
~104At.  The  Da-ci  values  that  were  averaged  over  the  period 
2  x  104At-4  x  104At  for  all  six  architectures  are  listed  in  Table  3. 

In  Fig.  5(b)  the  increase  in  Dci-ci  with  respect  to  architecture  I 
(A4[AC]-6),  is  plotted  against  DtopoU-Dtopol.b  where  Dt0poU  is 
defined  as  the  smallest  possible  topological  distance  between  two 
nearest  C  beads  (in  units  of  number  of  bonds  or  springs,  see 
Fig.  1(b)),  and  Dt0poi,2  is  the  smallest  possible  distance  between  a  C 
bead  and  its  2nd  nearest  C  bead.  ADci-ci  increases  non-linearly 
with  Dtopoi,2-Dtopoi,i-  These  results  are  in  line  with  studies  [39  in 
which  a  quadratic  dependency  was  found  for  in  total  18  Ax[AC]Ay 
[AC]  architectures  of  several  IEC  with  various  values  of  x  and  y.  In 
Ref.  39  the  bead  volumes  were  half  the  size  of  those  in  this  work 


The  pore  connectivity  can  be  probed  by  MC  trajectory  calcula¬ 
tions  by  restricting  the  tracer  particle  movement  to  the  pore  phase 
[34].  By  following  the  random  movement  of  the  particles,  tracer 
diffusion  coefficients  can  be  determined.  The  particle  movement 
takes  place  on  a  cubic  grid  of  size  2403  (1.4  x  107  nodes).  On  this 
grid  the  W  pore  network  is  constructed  such  that  nodes  for  which 
the  nearest  bead  is  a  W  bead  belong  to  the  pore  phase.  After  con¬ 
struction  of  the  pore  network  2000  particles  (NtraCer  =  2000)  are 
initially  put  at  randomly  selected  nodes  that  belong  to  the  pore 
phase.  During  every  Monte  Carlo  step  (MCS),  a  jump  trial  towards  a 
nearest  node  in  one  of  the  six  orthogonal  directions  is  randomly 
selected  for  each  particle  independently.  A  trial  is  successful  when 
the  aimed  site  belongs  to  the  W  (or  pore)  phase.  Diffusion  constants 
are  obtained  from  the  MSD  of  NtraCer  trajectories  (Eq.  (10)): 


(10) 


MSD  curves  obtained  for  morphologies  generated  at 
t  =  4  x  104At  are  displayed  in  Fig.  6(a).  Also  the  pure  water  case,  for 
which  every  jump  trial  is  accepted  and  the  slope  is  per  definition 
equal  to  1,  is  displayed  in  Fig.  6(a).  Diffusion  constants  within  the 
membranes,  DMg  can  be  expressed  relative  to  that  of  pure  water 
according  to  Eq.  (11)  [34,38-41  . 


A(MSD) 
MC  A  (MCS) 


(11) 


The  A(MSD)/A(MCS)  slopes  were  determined  by  linear  fitting 
over  the  time  window  3  x  105-6  x  105  MCS. 

Fig.  6(b)  displays  Dmc  values  against  DPD  time.  No  persistent 
increasing  or  decreasing  trend  of  DMc  with  DPD  time  is  observed 
beyond  2  x  104At.  This  indicates  that  equilibration  has  been 
reached,  and  is  in  accordance  with  conclusions  derived  from  g(r), 
where  Dci-ci  was  found  to  stabilize  beyond  ~104At  (Fig.  5(a)). 
Despite  the  apparent  convergence  of  Dmc  and  g(r),  Dmc  values  are 


4000  At 


8000  At 


12000  At 


16000  At 


20000  At 


24000  At 


28000  At 


32000  At 


Fig.  3.  Example  of  time  evolution  of  the  pore  network  for  A8[AC][AC]-3  architecture.  C  beads  are  colored  yellow  (light  grey)  and  W  beads  are  blue  (dark).  (For  interpretation  of  the 
references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 
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Fig.  4.  Time  evolution  ofW  bead  pair  correlation  functions  obtained  for  membranes  of  type  I:  A4[AC]-6  (a),  II:  A7[AC][Al[AC]-3  (b),  III:  A8[AC][AC]-3  (c)  and  IV:  A8[A2C2]-3  (d).  For 
increasing  time  the  curves  are  shifted  upwards  with  a  constant  along  the  vertical  axes. 


fluctuating  which  indicates  that  the  morphologies  are  continuously 
changing,  and  connections  between  pores  might  become  better 
(resulting  in  enhanced  DMc)  or  worse.  Diffusion  constants  through 
the  frozen  morphologies  averaged  over  the  time  window 
2  x  104Adt-4  x  104At  are  listed  in  Table  3  and  are  plotted  against 
Dci-ci  in  Fig.  7.  For  polymers  with  one  C  bead  (type  I,  II,  III)  within 
each  side  chain,  Dmc  increases  with  Dci-ci-  For  type  IV  polymer, 
which  contains  2C  beads  per  side  chain  this  is  also  the  case  with 
Dci-ci  being  substantially  larger. 


For  polymers  of  similar  length  with  Nbackbone  =  24,  diffusion 
turns  out  to  be  highest  for  architecture  III.  For  the  longer  IHb  and 
IVb  (Nbackbone  =  48)  architectures  water  diffusion  is  lower  than  for 
Ilia  and  IVa,  respectively  (Fig.  7).  Dci-ci  is  also  lower  for  the  longer 
architectures.  MC  diffusion  constants  were  obtained  from  frozen 
lattices,  therefore  differences  in  diffusion  are  caused  by  differences 
in  pore  morphologies.  Apparently  the  longer  architectures  tend  to 
be  less  capable  to  provide  good  connected  large  pores,  because 
more  C  beads  are  restricted  to  be  attached  to  the  same  polymer. 


♦  I:  A4[AC]-6 

□  II:  A7[AC]A1  [AC >3 
▲  111(a):  A8[AC][AC]-3 
a  1 11(b):  A8[AC][AC]-6 

•  IV(a):  A8[A2C2]-3 
o  IV(b):  A8[A2C2}-6 


0  10  20  30  40 

DPD  time  (lOOOAt) 


Fig.  5.  (a)  Dci-ci  plotted  against  DPD  time  for  all  six  architectures,  (b)  Increase  in  Dci-ci  (ADCi-ci  =  Dci-ci-Dci-ci(A4[AC]-6))  plotted  against  difference  in  topological  distance 
between  2nd  nearest  (Dt0p0i  2)  and  nearest  (Dt0poi,i)  C  beads  within  the  architectures.  Values  of  ADq-ci  are  averages  over  time  window  2  x  104At-4  x  104At.  The  drawn  curve  is  a 
quadratic  function  ( y  =  constant  xx2). 
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Table  3 

Average  distance  between  water  clusters  DCi-ci  and  diffusion  constants  DM c  derived 
from  MC  tracer  simulations. 


Architecture 

Da-ci 

Dm  c 

I:  A4[AC]-6 

4.25  nm 

0.048 

II:  A7[AC]Al[AC]-3 

4.62  nm 

0.085 

Ilia:  A8[AC][AC]-3 

5.02  nm 

0.131 

Illb:  A8[AC][AC]-6 

4.75  nm 

0.108 

IVa:  A8[A2C2]-3 

6.19  nm 

0.103 

IVb:  A8[A2C2]-6 

6.08  nm 

0.089 

I  conclude  with  the  following  remarks:  (i)  For  all  architectures 
Dcm  increases  with  Dt0poi,2-  L>topoi,i  (Fig.  5(b));  (ii)  For  both  [AC]  and 
[A2C2]  type  architecture  DMc  increases  with  DCi-ci  (Fig.  7);  (iii)  For 
architectures  (I,  II,  Ilia  and  IVa)  which  are  of  the  same  length  with 
^backbone  =  24  diffusion  is  highest  for  type  III  architecture.  Among 
the  architectures  for  which  Nbackbone  =  48  (Illb  and  IVb)  also  type  III 
reveals  the  highest  diffusion. 

3.3.  Largest  connected  water  cluster 

Whether  variations  in  Dmc  are  due  to  differences  in  the  amount 
of  water  being  located  within  isolated  clusters,  or  whether  diffusion 
is  to  various  extent  limited  by  the  presence  of  narrow  connections 
within  the  percolating  water  cluster  was  verified  by  calculating  the 
cluster  size  distribution  (CSD)  of  the  W  beads.  Here,  a  water  cluster 
is  defined  to  be  of  size  Nd  when  it  contains  Nd  W  beads  and  the 
distance  of  each  W  bead  within  that  cluster  to  at  least  one  of  the 
remaining  Nd- 1  beads  within  that  cluster  is  less  than  Lbond.  Due  to 
the  repulsive  interactions  between  DPD  beads  (Eq.  (2))  W  beads  are 
not  randomly  distributed  within  the  water  phase.  For  very  small 
Lbond  values  W  beads  can  thus  not  be  contained  within  large  clus¬ 
ters  and  each  W  bead  will  be  isolated.  But  for  Lbond  values  near  or 
well  above  the  first  peak  in  the  pair  correlation  function  (Fig.  4) 
water  clustering  occurs. 

Since  the  first  peak  in  the  pair  correlation  function  is  located 
near  0.85rc  W  clustering  is  expected  for  Lbond  near  or  above  0.85 rc. 
This  was  verified  by  calculating  the  CSD  as  function  of  Lbond  from 
the  positions  of  W  beads  for  an  equilibrated  DPD  system  that 
contained  only  W  beads.  The  size  of  the  largest  W  cluster  was  found 
to  increase  sharply  beyond  Lbond  =  0.8,  and  for  Lbond  =  0.85  more 
than  99.9%  of  the  total  number  of  W  beads  was  contained  in  one 
cluster.  For  all  morphologies  the  largest  W  cluster  was  obtained 
from  the  CSD  calculated  for  Lbond  =  0.85  and  Lbond  =  0.9.  The  results 
for  Lbond  =  0.85  are  plotted  against  DPD  time  in  Fig.  8.  Beyond  104At 
it  is  found  that  with  exception  of  architecture  I  (A4[AC]-6)  a  vast 
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Fig.  7.  Tracer  diffusion  coefficients  DMC  vs  inter  cluster  distance  Dq-ci  obtained  for  the 
architectures  listed  in  Table  l.The  backbone  lengths  are  given  in  parenthesis. 
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majority  of  the  total  amount  of  W  beads  is  located  within  a  single 
cluster.  The  largest  cluster  for  architecture  I  contains  between  20% 
and  40%  of  the  total  amount  of  W  beads  for  Lbond  =  0.85  (Fig.  8).  For 
Lbond  =  0.9  between  70%  and  90%  of  the  total  number  of  W  beads 
were  found  to  join  one  cluster  for  t  >  4  x  103At.  The  variations  in 
Dmc  observed  for  architectures  II,  III  and  IV  (Table  1 )  are  thus  not 
due  to  differences  in  percolating  cluster  size  but  due  to  more  salient 
differences  within  the  internal  structure  of  the  percolating  network 
(bottlenecks,  dead  ends  etc.). 


3.4.  DPD  W  bead  diffusion  through  dynamic  pore  networks 


Diffusion  coefficients  can  in  principle  also  be  calculated  from  the 
ensemble  average,  D^em^erffe t  of  the  MSD  of  W  beads  during  the 
DPD  simulations.  The  molecular  water  diffusion  coefficient  within 
the  membrane  D^2e^brane  is  then  obtained  by  comparing  £w^eadne 
with  the  W  bead  diffusion  coefficient  F)^irebe^dter  f°r  a  system  solely 
composed  of  W  beads  according  to  Eq.  (12a),  withD{£gwater  the 
pure  water  diffusion  coefficient  of  2.3  x  10-5  cm2  s-1. 


nMembrane 

r*Membrane  _  ^W-bead  _  rjPure  water 
H2O  —  jnPure  water  H20 

^W-Bead 


(12a) 


Ddpd 


nMembrane 

^W-bead 

nPure  water 
^W-Bead 


(12b) 


Instead  of  deriving  diffusion  coefficients  from  the  DPD  W  bead 
movements  the  MC  approach  was  applied  in  Section  3.2.  One  of  the 


♦  I :  A4[AC]-6  □  II:  A7[AC]A1  [ACT3 

a  I II (a):  A8[AC][AC]-3  a  I II (b):  A8[AC][AC]-6 

•  IV(a):  A8[A2C2]-3  o  IV(b):  A8[A2C2]-6 


MC  time  (105  MCS) 


— I:  A4[AC]-6 
-  -Q-  -  II:  A7[AC]A1  [ACP3 
— a —  111(a):  A8[AC][AC]-3 
A  1 11(b):  A8[AC][AC]-6 

•  IV(a):  A8[A2C2]-3 


DPD  time  (lOOOAt) 


Fig.  6.  (a)  Examples  of  MSD  curves  within  the  DPD  pore  networks  (t  =  4  x  104At)  for  architectures  I,  II,  Ilia,  Illb,  IVa  and  IVb.  The  pure  water  case  (0W  =  1.0),  for  which  the  slope  is 
equal  to  1,  is  included  (crosses),  (b)  DMc  plotted  against  DPD  time. 
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♦  I:  A4[AC]-6 

□  II:  A7[AC]A1  [AC]-3 
a  Ilia:  A8[AC][AC]-3 
Alllb:  A8[AC][AC]-6 

•  IVa:  A8[A2C2]-3 
olvb:  A8[A2C2]-6 


_ I 

0  10  20  30  40 

DPD  time  (1 03  At) 

Fig.  8.  Time  evolution  of  size  of  largest  water  cluster  for  all  six  architectures, 
i-bond  =  0.85. 

reasons  was  that  DPD  allows  bond  crossings  that  lift  up  entangle¬ 
ments.  This  causes  DPD  polymers  to  diffuse  faster  than  real  poly¬ 
mers,  resulting  in  fast  evolvement  to  equilibrium  like 
morphologies.  As  a  result,  the  DPD  pore  morphology  may  change 
while  the  water  beads  are  diffusing  through  it.  For  this  reason, 
frozen  (static)  DPD  morphologies  were  taken  as  input  for  the  MC 
trajectory  calculations  and  the  pore  networks  were  thus  fixed  in 
time. 

Since  the  pores  can  rearrange  continuously  the  relative  W  bead 
diffusivities,  Dd pd  (defined  by  Eq.  (12b)),  expressed  relative  to  pure 
water  diffusion  are  expected  to  be  higher  than  the  DMc  values.  The 
DDpd  values  are  compared  with  Dmc  in  Fig.  9(a)  (filled  symbols).  Dmc 
and  DDpd  are  highest  for  type  III  architectures  and  lowest  for  ar¬ 
chitecture  I.  As  expected,  the  Ddpd  diffusivities  differ  from  the  MC 
derived  diffusivities  with  Ddpd  being  between  a  factor  2.5  and  4.3 
higher  than  the  Dmc  values. 

C  beads  are  the  end  moieties  of  the  pendant  side  chains 
(Fig.  1(b))  and  are  located  near  the  W  pore  network  (Fig.  2).  Their 
pendant  motion  (and  A  bead  motions  as  well)  result  in  a  continuous 
rearrangement  of  the  W  pore  network  (Fig.  3).  Now  naively  assume 
that  the  A  beads  are  fixed  in  space  for  t  >  20000At.  The  W  beads  can 
then  have  access  to  at  most  the  whole  hydrophilic  phase  which 
constitutes  a  volume  fraction  of  0.30  (</>w+c  =  0w+^c  =  0.3). 
Additional  MC  trajectory  calculations  were  repeated  within  the 
mapped  DPD  morphologies  sampled  at  t  >  20000At  in  which  tracer 
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particles  can  diffuse  freely  through  the  mapped  hydrophilic 
(W  +  C)  networks.  Dd pd  is  plotted  against  thus  calculated  MC 
diffusion  constants,  Dmc(W  +  C)  in  Fig.  9(b).  Dd pd  is  systematically 
higher  than  DMc(W  +  C).  W  beads  are  thus  more  mobile  than  ex¬ 
pected  for  a  fixed  C  +  W  pore  network.  Apparently  the  mobility  of 
the  A  beads  also  causes  increase  of  W  bead  diffusion  due  to  the 
spatiotemporal  evolvement  of  the  pore  networks. 

During  a  DPD  simulation  the  C  +  W  pore  network  can  be  made 
nearly  static  by  increase  of  the  polymer  A  bead  masses  which  is 
expected  to  result  in  a  lower  W  bead  mobility.  This  was  verified  by 
rerunning  the  DPD  simulations  starting  from  the  configurations 
shown  in  Fig.  2  (t  =  2  x  104At).  Dd pd  obtained  for  ma-bead  =  1000 
and  mc-bead  =  m^-bead  =  1  are  given  in  Fig.  9(a).  Diffusivities  of  the 
DPD  W  beads  turn  out  to  be  still  a  factor  ~1.3 — 1.7  larger  than  the 
DMc  averages.  A  further  decrease  of  Dd pd  is  obtained  by  also 
increasing  the  C  bead  mass.  For  uia- bead  =  mc-bead  =  1000  the  Dd pd 
values  approach  those  obtained  from  Dmc  for  architecture  I  but  for 
the  other  architectures  are  less  than  DMc  (Fig-  9(a)).  Apparently  the 
decreased  mobility  of  the  C  beads  strongly  influences  W  bead 
transport  through  the  nearly  static  pore  networks.  It  is  well  possible 
that  the  low  C  bead  motion  for  mc-bead  =  1000  causes  W  beads  to  be 
caged  in  and  trapped  for  longer  periods  resulting  in  a  significant 
decrease  of  Dd pd- 


4.  Discussion 

The  polymer  lengths  for  architectures  I,  II,  Ilia  and  IVa  were  the 
same  (Nbackbone  =  24).  The  calculated  tracer  diffusion  constant,  Dmc, 
for  water  diffusing  through  the  frozen  W  network  is  largest  for 
architecture  Ilia  (Fig.  7).  When  tracer  particles  are  also  allowed  to 
move  through  the  W  +  C  network  diffusion  increases  but  the  trends 
are  not  affected  with  architecture  Ilia  revealing  the  highest  diffu¬ 
sion  constant,  Dmc(W  +  C)  (Fig.  9(b)).  For  each  architecture 
Dmc(W  +  C)  is  only  slightly  less  than  the  W  bead  diffusion  constant 
Ddpd  (Fig.  9(b)),  which  suggested  that  W  beads  diffusing  through 
the  “spatiotemporally  evolving  membrane”  can  pass  through  the 
region  previously  occupied  by  a  C  bead.  In  the  DPD  simulations,  the 
polymer  motions  thus  promote  W  bead  transport.  This  was  verified 
by  decreasing  the  A  or  (A  and  C)  bead  mobility  by  increasing  their 
masses.  This  caused  a  significant  lower  W  bead  diffusion  which 
approach  the  MC  values  for  diffusion  through  the  frozen  W  net¬ 
works  (Fig.  9(a)).For  real  membranes  a  reduction  of  polymer  mo¬ 
tion  occurs  by  temperature  decrease,  which  is  expected  to  reduce 
polymer  assisted  water  diffusion.  Hypothetically,  the  Dmc  diffusion 
constants  obtained  for  tracer  diffusion  through  the  static  DPD  pore 
networks  are  then  lower  bounds  for  low  temperature  conditions 
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Fig.  9.  (a)  DPD  W  bead  diffusion  constants  vs  DMc  within  the  static  DPD  morphologies.  Bead  masses  are  equal  to  1  (filled  diamonds):  mA_bead  =  103  (open  triangles):  mA_bead  =  rnC- 
bead  =  103  (open  circles),  (b)  Ddpd  vs  DMc(W  +  C).  All  bead  masses  are  equal  to  1. 
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(when  morphology  is  not  affected  and  locally  water  diffusion  re¬ 
sembles  that  of  pure  water  at  that  lower  temperature). 

Comparison  between  the  larger  polymers  Illb  and  IVb  with 
^backbone  =  48,  also  reveals  for  type  III  polymer  the  highest  water 
diffusion  constants,  Dmc  (Figs.  9(a)  and  7),  Dmc(W  +  C)  (Fig.  9(b)) 
and  Ddpd  (Fig.  9(a)).  For  the  longer  Illb  and  IVb  sequences  the  W 
bead  diffusion  constant  Dd pd  (Fig.  9)  is  lower  than  for  Ilia  and  IVa, 
respectively.  Since  also  within  the  static  morphologies  Dmc  (Fig.  7) 
and  DMc(W  +  C)  (Fig.  9)  are  lower  for  the  longer  architectures,  the 
difference  is  not  attributed  to  a  different  extent  to  which  polymer 
motions  assist  in  W  bead  diffusion.  The  decrease  in  W  bead  diffu¬ 
sion  is  thus  due  to  different  pore  morphologies  with  lower  DCi-ci 
values  for  the  longer  (Illb  and  IVb)  architectures  (see  Fig.  7).  Within 
real  polymer  membranes  intrinsic  backbone  stiffness,  due  to 
bending  constraints  along  the  polymer  backbone,  may  also  affect 
the  formation  of  good  connected  pores  required  for  fast  water 
diffusion.  In  the  present  simulations  bending  stiffness  was  not 
included. 

Water  diffusion  derived  from  MC  simulations  through  the 
frozen  (W  or  W  +  C)  networks  shows  the  same  trends  with  respect 
to  polymer  architecture  as  derived  from  the  DPD  W  bead  motions 
through  the  evolving  pore  network  (see  Fig.  9).  Trends  are  also  not 
affected  by  decreasing  the  polymer  mobility  by  selecting  higher 
polymer  bead  masses.  This  validates  to  derive  correlations  between 
polymer  architecture  and  water  transport  in  the  membrane.  This  is 
because  water  only  diffuses  through  the  hydrophilic  (W  or  W  +  C) 
phase  and  in  both  MC  and  DPD  modeling  the  local  water  diffusion 
constant  in  the  water  phase  is  the  same  as  that  of  pure  water.  In 
modeling  water  transport,  the  use  of  frozen  DPD  morphologies  as 
input  for  MC  tracer  diffusion  calculations  to  predict  relations  be¬ 
tween  polymer  architecture  and  water  transport  is  therefore  an 
appropriate  approach. 

In  modeling  gas  transport  more  care  should  be  taken  into  ac¬ 
count.  Within  membranes  the  gas  concentration  and  mobility  de¬ 
pends  on  the  phase  in  which  the  gas  species  resides.  Gas 
permeation  cannot  be  simulated  by  conventional  DPD  since  it  re¬ 
quires  the  introduction  of  gas  beads  that  reproduce  gas  solubility 
and  gas  diffusion  within  polymer  and  water  phase.  DPD  in  combi¬ 
nation  with  (kinetic)  MC,  however,  can  reveal  trends  that  relate  the 
EW  and  water  content  dependence  of  O2  permeability  through 
Nation  [37,39].  These  studies  rely  on  the  assumption  that  the  gas 
solubility  and  diffusion  in  both  phases  are  the  same  as  in  the  two 
pure  component  phases.  O2  permeation  in  Nation  has  also  been 
proposed  to  be  promoted  by  polymer  flexibility  53]  and  to  occur 
mainly  in  the  intermediate  region  located  in  the  polymer  matrix 
near  the  ionic  clusters  [53,54].  This  cannot  be  modeled  by  DPD 
since  it  would  require  the  reproduction  of  the  gas  solubility  and 
diffusion  within  three  distinct  (water,  polymer,  and  intermediate) 
phases.  Also  KMC  calculation  of  gas  permeation  through  frozen 
lattices  obtained  from  DPD  requires  these  six  parameters  as  input. 
Only  when  these  are  known  (from  experiment  or  MD  simulations) 
DPD  combined  with  KMC  might  be  value  to  model  gas  permeation 
through  meso-scale  pore  networks. 

The  interactions  between  hydrophobic  (A)  and  hydrophilic  (C, 
W)  beads  were  chosen  close  to  those  used  in  previous  modeling  on 
Nation.  Therefore  it  may  be  expected  that  membranes  composed  of 
PFSA-like  polymers  of  type  III  (Fig.  1(b))  result  in  better  pore  con¬ 
nectivity  and  higher  solvent  diffusion  constants  and  proton  con¬ 
ductivity  as  well.  As  an  example,  two  additional  simulations  were 
performed  using  the  x-parameters  proposed  for  Nation  in  Ref.  42 
(given  in  Section  2.2).  According  to  Eq.  (7)  the  DPD  repulsions  are: 
Qa’B’  =  104.1 ;  Qa’c  =  116.9,  Qa’w  =  128.1;  ub’c  =  109.7;  ob’w  =  124.4; 
acw  =  92.4.  Using  these  repulsions  the  obtained  values  for  Dmc  at 
t  =  20000A t  for  A4'[B'C']-6  and  AS'IB'C'HB'C'J-S  at  </>w  =  0.16  are 
0.014  and  0.122,  respectively.  For  DMc(W  +  C)  respectively  0.064 


(A4'[B'C']-6)  and  0.246  (AS'lB'C'nB'C'1-3)  are  obtained.  Diffusion 
within  AS'IB'C'HB'C']  membranes  are  thus  systematically  higher 
than  within  A4'[B/C/]  membranes. 

When  the  DMc(W  +  C)  values  are  converted  to  physical  diffusion 
constants  Dphys(H20)  then  1.5  x  1CT6  cm2  s-1  and 
5.5  x  10~6  cm2  s^1  are  obtained  for  the  A4'[B'C']-6  and  AS'lB'C'] 
[B'C'J-3  architectures,  respectively.  For  comparison  the  experi¬ 
mental  value  for  NafionllOO  at  0W  =  0.16  derived  from  Fig.  4(a)  in 
Ref.  11  is  ~2  x  10-6  cm2  s-1,  which  is  in  the  same  range  as  derived 
here. 

The  above  example  illustrates  that  enhanced  water  diffusion 
predicted  for  two  side  chains  branching  off  from  backbones  does 
not  only  occur  for  the  repulsions  used  in  Table  2  but  for  a  range  DPD 
repulsions.  PFSA  polymers  in  a  PEFC  may  contain  several  hundred 
side  chains  [55],  and  their  distribution  along  the  backbone  is  not 
precisely  known.  Also  the  dependency  of  water  diffusion  on  poly¬ 
mer  length  is  not  known.  An  exact  reproduction  of  experimental 
diffusion  constants  can  thus  not  be  expected. 

In  total  120  pore  morphologies  that  were  obtained  by  DPD  for 
six  architectures  were  subjected  to  MC  tracer  diffusion  simulations 
on  a  lattice.  These  MC  simulations  were  robust,  neglecting  any 
dynamical  coupling  between  water  and  polymer,  and  no  correla¬ 
tions  between  the  water  molecules  were  assumed.  The  only 
assumption  was  that  the  local  molecular  water  mobility  within  the 
pores  is  the  same  as  that  within  pure  water.  Despite  this  over¬ 
simplification,  quasi  elastic  neutron  scattering  studies  [56,57]  on 
hydrated  Nation  membranes  indeed  concluded  that  within  the  nm 
scale  pores  the  water  diffusion  constant  is  close  to  that  of  pure 
water.  For  instance,  Pivovar  and  Pivovar  [57]  measured  at  a  water 
content  of  ~14  vol%  a  local  water  mobility  of  1.6  x  10  5  cm2  s-1, 
which  is  only  30%  below  the  pure  water  diffusion  constant  of 
2.3  x  10-5  cm2  s-1.  Since  the  obtained  MC  tracer  diffusion  constants 
within  the  (W  or  W  +  C)  pores  correlate  strongly  with  the  DPD  W 
bead  diffusion  constant,  it  is  obvious  that  it  is  the  pore  topology 
that  mostly  affects  the  diffusion  of  water  beads  (molecules). 

Flexible  polymer  backbones  and  side  chains  were  assumed  (no 
bending  constraints).  This  allowed  the  relaxation  towards  mor¬ 
phologies  with  optimized  compatible  (W-W,  C-W,  A-A)  contacts 
and  minimized  number  of  incompatible  hydrophobic-hydrophilic 
(A-W,  A-C)  contacts.  Backbone  chain  lengths  as  reported  for 
realistic  PFSA  polymers  [55]  are  at  least  one  order  of  magnitude 
larger  than  here.  Further  improvement  of  diffusion  might  be  to  see 
if  there  is  for  each  architecture  type  an  optimal  polymer  length 
which  results  in  highest  water  diffusion.  For  the  sake  of  a  fair 
comparison,  the  polymer  backbone  lengths  for  architectures  I,  II, 
Ilia,  and  IVa  were  intentionally  kept  fixed  at  Nbackbone  =  24.  But  a 
dependency  on  polymer  length  was  found  when  comparing  ar¬ 
chitecture  Ilia,  IVa  vs.  Illb,  IVb  (Nbackbone  =  48).  The  effect  of  intrinsic 
chain  stiffness,  not  considered  here,  may  affect  the  pore 
morphology  and  water  diffusion.  For  fixed  chain  stiffness  an  opti¬ 
mum  in  chain  length  might  be  present  to  facilitate  water  (cq. 
proton)  diffusion.  Also  adjusting  x  parameters  might  optimize 
diffusion.  For  block  polymers  also  a  ^-parameter  dependency  on 
diffusion  was  found  [35].  Fine  tuning  by  searching  for  optimized 
polymer  lengths,  ^-parameters,  or  bending  stiffness  was  not  the 
intention  of  this  work,  but  worthwhile  to  look  for  in  future  studies. 

For  realistic  membranes  that  are  applied  in  a  PEFC,  some  rigidity 
is  necessary.  The  rigidity  is  partly  obtained  by  the  presence  of  micro 
crystals  [11].  DPD  does  not  permit  to  study  polymer  crystallization. 
Therefore  the  extent  to  which  crystallization  impairs  the  formation 
of  the  optimized  morphologies  could  not  be  studied  here.  It  is  well 
possible  that  the  amount  of  crystallization  is  higher  for  architec¬ 
tures  with  side  chains  separated  far  apart  along  the  polymer 
backbone.  In  the  present  study  these  are  the  architectures  of  type  III 
(A8[AC][AC])  and  IV  (A8[A2C2])  for  which  also  the  largest  Dci-ci 
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were  obtained  (Fig.  7).  Interestingly,  the  larger  hydrophobic  space 
between  the  clusters  for  these  architectures  as  compared  to  A4[AC] 
architectures  might  favor  the  formation  of  small  polymeric  crystals 
necessary  for  the  required  rigidity. 

Electrostatic  interactions  were  not  explicitly  taken  into  account 
in  the  DPD  simulations.  Therefore  the  phase  separated  morphol¬ 
ogies  in  this  work  are  a  consequence  of  both  bead  connectivity 
(chain  architecture)  and  the  choice  of  %  parameters.  It  would  be 
interesting  to  see  to  what  extent  the  results  are  affected  when 
electrostatic  interactions  are  explicitly  taken  into  account.  From 
earlier  DPD  studies  [39,40]  it  was  found  that  C  beads  which  are 
topologically  close  together  within  a  certain  architecture  (low  Dt0_ 
pol  l  value)  tend  to  join  the  same  pore  or  water  cluster.  For  C  beads 
that  are  close  to  each  other,  the  electrostatic  interactions  between 
them  are  also  expected  to  affect  the  distance  between  them.  For 
low  hydration  levels  the  C  beads  have  a  neutral  character.  With 
increasing  hydration  proton  dissociation  occurs,  which  leaves  the 
protons  in  the  water  clusters.  This  may  result  in  repulsion  between 
nearby  C  beads  with  increase  in  hydration  level.  These  issues  were 
beyond  the  scope  of  this  work  and  should  be  challenged  by  ab  initio 
or  MD  methods.  Noteworthy  is  that  a  total  reflectance  infrared 
spectroscopy  study  [58]  revealed  that  for  Nation  the  sulfonic  site 
is  protonated  at  dry  and  dissociated  in  hydrated  Nation. 
Combined  with  ab  initio  calculations  Webber  et  al.  [58]  found  that 
sulfonic  acid  sites  are  deprotonated  above  a  hydration  level  of  A  =  4, 
in  accordance  with  other  ab  initio  work  on  Nation  or  Dow  mem¬ 
branes  [59  . 

One  of  the  main  conclusions  is  that  the  simulations  presented  in 
this  work  predict  that  ionomers  composed  of  A8[AC][AC]  and  A'8 
[B'C'HB'C']  repeat  units  give  rise  to  faster  water  diffusion,  or  better 
connected  pore  networks,  as  compared  to  those  composed  of  A4 
[AC]  and  A'4[B'C']  repeat  units,  respectively.  To  my  knowledge,  at 
present  it  is  difficult  to  validate  this  prediction  experimentally  since 
precise  control  of  branching  distribution  is  not  feasible  during 
polymer  synthesis. 

Finally,  it  would  be  worthwhile  to  perform  a  fine  graining  of 
the  present  work  for  system  size  of  O(104  nm3)  by  MD  simula¬ 
tions  for  the  architectures  in  Fig.  1(b).  As  mentioned  in  the 
Introduction  the  evolvement  from  pre-assumed  morphologies 
[31  ]  can  then  be  studied  for  such  system  size.  After  fine  graining 
the  equilibrium-like  DPD  morphologies,  they  can  then  serve  as 
starting  morphologies  to  study  their  stability  and  evolvement. 
Also  the  evolvement  from  random  configurations  toward  equi¬ 
librium  structures  might  be  studied  by  MD.  By  repeatedly 
creating  periodic  images  that  are  again  equilibrated  large  scale 
structures  can  then  be  obtained  [60  .  Support  for  the  findings 
here  can  then  be  obtained  when  a  comparison  is  made  between 
the  molecular  water  diffusion  constants  and  that  of  W  beads 
through  the  DPD  morphologies  for  the  architectures  in  Fig.  1(b). 
Validation  requires  that  similar  trends  with  respect  to  chain  ar¬ 
chitecture  are  obtained  for  both  fine  grained  (MD)  and  coarse 
grained  (DPD)  simulation. 

5.  Conclusion 

Pore  morphologies  within  meso  phase  separated  hydrated 
grafted  polymeric  membranes  were  simulated  by  DPD  in  which 
polymer  fragments  are  represented  by  connected  beads.  Water 
diffusion  through  the  water  containing  pores  was  modeled 
separately  by  means  of  MC  tracer  calculations  through  selected 
(static)  morphologies.  Four  types  of  membranes  were  considered 
that  are  of  equal  IEC  but  differ  in  the  way  side  chains  are  attached 
to  the  hydrophobic  backbones.  The  main  and  side  chains  were 
flexible  with  no  bending  stiffness  taken  into  account,  therefore 
the  phase  separated  pore  morphologies  are  solely  governed  by 


the  adapted  Flory-Huggins  %  parameters,  polymer  architecture, 
and  water  contents  (16  volume  %).  With  A  beads  being  hydro- 
phobic  and  C  beads  hydrophilic,  it  was  found  that  for  polymers 
composed  of  similar  length  of  24  backbone  beads  the  MC  derived 
diffusion  constants  through  static  (frozen)  pore  networks  in¬ 
crease  as  A4[AC]  <  A7[AC]A1[AC]  <  A8[A2C2]  <  A8[AC][AC]. 
Simulations  performed  for  the  A8[A2C2]  and  A8[AC][AC]  archi¬ 
tectures  of  two  times  longer  length  (48  backbone  beads), 
revealed  again  the  A8[AC][AC]  polymer  showing  the  highest 
water  diffusion  constant  of  both. 

The  results  are  complemented  by  calculation  of  diffusion  of  the 
water  beads  through  the  evolving  dynamic  DPD  pore  network.  The 
same  trend  as  derived  from  the  MC  calculations  was  obtained  with 
water  diffusion  being  systematically  higher  for  the  A8[AC][AC]  ar¬ 
chitecture.  This  trend  turned  not  be  affected  when  the  mobility  of 
the  polymer  A  and/or  C  beads  was  artificially  decreased  by  increase 
of  the  bead  masses.  These  new  insights  suggests  that  within  (new 
types  of)  membranes  of  equal  IEC  but  with  two  side  chains 
branching  off  from  a  single  backbone  site  allows  the  formation  of 
better  connected  hydrophilic  pores. 
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